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NON-IDEAL  GAS  EFFECTS  ON  SHOCK  WAVES  IN  WEAKLY  IONIZED  GASES 

ROBERT  RUBINSTEIN*  AND  AARON  H.  AUSLENDERt 


Abstract.  A  shock  wave  in  a  weakly  ionized  gas  can  be  preceeded  by  a  charge  separation  region  if  the 
Debye  length  is  larger  than  the  shock  width.  It  has  been  proposed  that  electrostatic  contributions  to  pressure 
in  the  charge  separation  region  can  increase  the  sound  speed  ahead  of  the  shock  well  above  the  sound  speed 
in  a  neutral  gas  at  the  same  temperature  and  therefore  increase  the  shock  propagation  speed.  This  proposal 
is  investigated  numerically  and  theoretically.  It  is  concluded  that  although  the  ion  gas  becomes  strongly 
non-ideal  in  the  charge  separation  region,  there  is  no  appreciable  effect  on  the  neutral  shock. 

Key  words,  plasma  drag  reduction,  weakly  ionized  gases 

Subject  classification.  Fluid  Mechanics 

1.  Introduction.  Ionized  fluid  systems  are  effectively  electrically  neutral  at  length  scales  larger  than 
the  Debye  length ,  the  scale  at  which  charge  screening  becomes  important  [1].  In  shock  waves  in  weakly  ionized 
gases,  the  great  mass  disparity  between  electrons  and  ions  promotes  the  formation  of  a  charge  separation 
region  near  the  neutral  shock  front.  The  size  of  this  region  is  of  the  order  of  the  Debye  length;  it  might 
modify  the  shock  dynamics  if  the  Debye  length  is  larger  than  the  shock  width  [2]. 

Experimental  observations  [3]  suggest  that  shock  propagation  in  weakly  ionized  gases  is  ‘anomalous’  in 
some  respects;  for  example,  shocks  appear  to  travel  faster  than  they  would  in  a  completely  neutral  gas  at 
the  same  temperature  and  pressure.  It  has  been  proposed  [4]  that  the  unscreened  electrostatic  forces  in 
the  charge  separation  region  could  cause  a  significant  increase  in  pressure  ahead  of  the  shock,  leading  to  the 
sound  speed  increase  thought  to  be  responsible  for  these  effects.  This  pressure  increase  was  at  first  estimated 
to  be  quite  large  [4];  however,  it  cannot  simply  be  added  to  the  neutral  pressure.  Only  a  weaker  conclusion 
is  possible,  namely  that  electrostatic  effects  cause  the  ion  and  electron  gases  to  become  strongly  non-ideal 
in  the  charge  separation  region. 

Any  mechanism  which  increases  shock  propagation  speeds  requires  an  energy  source  in  the  medium 
into  which  the  shock  propagates.  In  this  proposal,  the  energy  source  is  the  electrostatic  potential  energy  in 
the  charge  separation  region.  This  possibility  is  attractive,  because  the  existence  of  the  charge  separation 
region  requires  a  large  Debye  length,  hence  a  relatively  low  density  of  charged  particles.  If  the  ‘anomalous 
effects’  are  real  and  are  associated  with  the  existence  of  the  charge  separation  region,  then  the  effects  must 
vanish  both  in  the  limit  of  very  high  charged  particle  density,  corresponding  to  a  small  Debye  length,  and 
in  the  limit  of  very  low  density,  corresponding  to  a  large  Debye  length  but  no  charged  particles,  hence  to  an 
effectively  purely  neutral  shock.  This  would  explain  the  absence  of  ‘anomalous  effects’  on  re-entry  vehicles 
and  suggests  that  an  optimal  density  might  exist  at  which  the  anomalous  effects  are  most  pronounced. 

This  possibility  is  investigated  numerically  and  theoretically.  Increased  ion  pressure  alone  is  found  to 
have  no  effect  on  the  neutral  shock,  but  the  shock  can  be  accelerated  by  an  increased  neutral  pressure, 
even  if  the  increase  only  occurs  in  the  charge  separation  region.  However,  because  the  electrostatic  force 
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responsible  for  the  increased  ion  pressure  does  not  act  on  the  neutral  particles,  this  pressure  increase  cannot 
be  transferred  to  the  neutrals.  Consequently,  the  increased  ion  pressure  in  the  charge  separation  region 
ahead  of  the  shock  front  in  a  weakly  ionized  gas  will  not  significantly  enhance  the  neutral  shock  propagation 
speed. 


2.  Governing  Equations  for  a  Shock  in  a  Weakly  Ionized  Gas.  Following  previous  work  on 
shock  waves  in  weakly  ionized  gases  [2],  [5],  [6],  the  weakly  ionized  gas  will  be  treated  as  a  three- component 
system  consisting  of  neutrals,  ions,  and  electrons.  The  continuity  equations  are 


(2.1) 

(2.2) 

(2.3) 


dn  d(nu)  _ 
dt  dx 
drii  d{niUi)  _ 
dt  dx 

dne  d(neue)  _ 
dt  dx 


where  R  denotes  ionization  and  recombination.  The  momentum  equations  are 
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and  the  energy  equations  are 


(2.7) 

(2.8) 
(2.9) 
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If  each  component  were  an  ideal  gas,  then  the  equations  of  state  would  be 


(2.10) 

p  =  nkT 

(2.11) 

pi  zz  mkTi 

(2.12) 

Pe  —  riekTe 

Modified  equations  of  state  corresponding  to  non-ideal  gas  effects  will  be  considered  later.  The  system  of 
governing  equations  is  completed  by  the  one-dimensional  relation 


(2.13) 


dE 

dx 


—  {rii  -ne) 


between  the  electric  field  and  the  density  of  charged  particles;  this  equation  is  equivalent  to  Poisson’s  equation 
for  the  electrostatic  potential. 

In  Eqs.  (2.1)— (2.13),  variables  without  a  subscript  pertain  to  the  neutral  particles,  the  subscripts  i  and 
e  denote  ions  and  electrons  respectively,  n  is  the  number  of  particles,  u  the  velocity,  p  the  pressure,  T  the 
temperature,  E  the  electric  field,  and  m  the  mass  of  a  particle,  so  that  the  product  mn  is  the  mass  density. 
The  masses  of  ions  and  neutrals  have  been  assumed  equal.  Viscous  and  other  additional  terms  are  included 
in  the  final  V  terms  in  Eqs.  (2.4)— (2.9);  all  of  these  terms  are  given  explicitly  in  [6]. 
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The  experimental  conditions  under  which  anomalous  effects  have  been  reported  are  conditions  of  thermal 
non-equilibrium  with  Te  T*  =  T  ahead  of  the  shock.  Maintaining  this  condition  in  the  presence  of 
collisional  coupling  requires  a  source  term  in  the  electron  energy  equation  [7].  Behind  the  shock,  the  system 
relaxes  to  a  state  of  thermal  equilibrium  with  Te  =  —  T.  For  the  analogous  problem  of  relaxation  behind 

shocks  in  vibrationally  excited  media,  see  [8].  The  possible  role  of  thermal  nonequilibrium  ahead  of  the  shock 
on  the  decay  of  anomalous  properties  has  been  considered  earlier  by  us  [9]. 

The  most  important  property  of  these  equations  is  that  the  coupling  between  neutrals  and  charged 
species  is  through  the  collisional  terms  vin  and  ven  only.  Numerical  solution  of  these  equations  [5]  do  not 
reveal  any  anomalous  effects  whatsoever  because  of  the  weakness  of  these  couplings.  The  present  work  will 
investigate  the  possibility  raised  in  [4]  that  by  replacing  the  equations  of  state  Eqs.  (2.10)-(2.12)  by  modified 
equations  of  state  which  reflect  the  non-ideal  gas  behavior  in  the  charge  separation  region,  anomalous  shock 
propagation  speed  can  be  predicted. 

3.  Numerical  Simulations.  The  governing  equations  are  integrated  using  a  code  developed  by  Adam¬ 
ovich  et  ah  [6].  This  code  is  based  on  a  general-purpose  solver  for  second-order  systems  of  partial  differential 
equations:  pdecol  [10]  which  integrates  systems  of  the  form 


(3.1) 


-  Fft  U  —  — ) 
dt  F(*,u’  dx’  dx2’ 


Fourth-order  spatial  spline  interpolation  is  used  with  continuity  of  second-order  derivatives  required  at  each 
node.  Because  the  solution  method  requires  two  boundary  conditions  for  each  equation,  artificial  viscous 
terms  are  added  to  the  continuity  equations  Eqs.  (2.1)— (2.3).  The  gridpoints  are  chosen  nonuniformly  with 
tank  stretching  to  place  the  finest  spatial  resolution  inside  the  shock.  The  solution  domain  is  the  finite 
normalized  region  0  <  x  <  1;  boundary  conditions  at  the  endpoints  are  assigned  based  on  a  Mach  2  neutral 
shock  with  conditions  on  the  charged  species  set  for  representative  experimental  conditions  [6].  Ahead  of 
the  shock,  a  state  of  thermal  nonequilibrium  Te  >  7*  =  Te  exists,  which  relaxes  to  thermal  equilibrium  with 
Te  —  Ti=T  behind  the  shock. 

The  shock  motion  is  computed  over  a  predetermined  time  step,  then  the  spatial  coordinates  are  recalcu¬ 
lated  by  setting  the  shock  center,  identified  by  the  position  of  maximum  gradient  of  the  neutral  density,  to 
the  reference  value  x  =  0.5.  Consequently,  if  the  shock  begins  to  move  too  quickly,  the  calculation  becomes 
inaccurate  because  the  boundary  conditions  are  contaminated  by  the  quantities  near  the  shock. 

In  order  to  obtain  results  rapidly  for  a  large  number  of  test  cases,  the  code  was  run  at  low  spatial 
resolution,  with  200  space  points. 

A  preliminary  code  verification  was  preformed  by  computing  a  neutral  shock,  obtained  by  setting  the 
collisional  couplings  to  zero.  It  was  found  that  the  predicted  motion  is  consistent  with  the  postulated 
boundary  conditions  for  a  Mach  2  shock.  A  simple  test  of  the  inability  of  the  equations  of  motion  Eqs. 
(2.1)-(2.13)  to  capture  anomalous  effects  is  provided  by  restoring  the  collisional  couplings.  It  was  found 
that  the  effect  on  the  shock  path  is  negligible.  Computational  results  illustrating  this  point  will  be  discussed 
below. 


4.  Simulation  Results. 

4.1.  The  charge  separation  region.  The  density  profiles  across  the  shock  are  shown  in  Fig.  (4.1) 
for  the  charge  densities  a  =  1CT5, 10“6  and  10-8  in  the  left,  center,  and  right  graphs  respectively.  The  value 
a  —  10-6  is  believed  to  be  typical  of  experimental  conditions. 
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density 


Fig.  4.1.  neutral  (N),  ion  (I),  and  electron  (E)  densities  in  shock  with  a  —  (left)  10  5;  (center)  10  6;  (right)  10  8.  The 
size  of  the  charge  separation  region  varies  as  a""1/2.  The  conditions  in  the  center  graph  are  believed  typical  of  experimental 
conditions. 


The  calculations  show  a  small  but  significant  charge  separation  region  when  a  =  10-6.  In  the  units  used 
in  the  calculation,  the  Debye  length  lD  =  0.0025,  which  is  consistent  with  the  width  «  0.005  =  21 D  of  the 
charge  separation  region  in  Fig.  (4.1).  The  charge  separation  region  disappears  completely  when  a  =  10~5, 
but  becomes  much  larger  when  a  =  10“8.  Again,  since  ld  ~  a"1/2,  ld  =  0.025  in  these  units,  and  the  region 
of  unbalanced  positive  charge  has  width  «  21  p. 

Note  that  the  densities  are  all  normalized  in  these  figures  to  have  the  same  dimensionless  values  ahead 
of  and  behind  the  shock;  thus,  although  the  charge  separation  region  is  very  large  when  a  =  10-8,  it  must 
be  understood  that  the  density  of  charged  particles  is  extremely  small  in  this  case. 

4.2.  The  effect  of  ion  pressure.  We  next  investigate  whether  an  increased  ion  pressure  alone  can 
change  the  shock  path.  In  the  a  =  1CT6  case,  multiplying  the  ion  pressure  gradient  by  factors  as  large  as  50 
has  no  effect  on  the  shock  path.  The  explanation  appears  to  lie  in  the  very  weak  collisional  coupling  and  the 
low  ion  density.  But  increasing  the  neutral  pressure  gradient  by  as  little  as  25%  causes  a  significant  increase 
in  shock  propagation  speed.  These  results  are  illustrated  in  Fig.  (4.2)  which  compares  the  shock  paths  for  a 
neutral  shock,  the  shock  with  increased  ion  pressure  gradient,  and  the  shock  with  increased  neutral  pressure 
gradient.  The  plotted  paths  are  the  changes  about  the  path  of  a  neutral  Mach  2  shock.  This  figure  illustrates 
the  point  made  above,  that  coupling  to  the  charged  species  has  no  effect  on  the  shock  path,  even  if  the  ion 
pressure  gradient  is  greatly  enhanced.  Thus,  shock  path  modification  requires  an  increase  in  the  neutral 
pressure. 

4.3.  The  effect  of  neutral  pressure.  To  test  the  effect  of  modifying  the  neutral  equation  of  state  in 
the  charge  separation  region,  the  neutral  equation  of  state  Eq.  (2.10)  was  replaced  by 

(4.1)  p  =  nkT{l  +  (7max[ni  -ne,0]} 

where  the  term  proportional  to  C  represents  an  enhancement  of  the  neutral  pressure  due  to  charge  imbalance. 
The  specific  mechanism  for  this  enhancement  will  be  considered  later;  the  goal  of  this  section  is  to  evaluate 
the  possibility  of  shock  propagation  speed  enhancement  through  a  modified  neutral  equation  of  state.  The 
constant  C  was  given  the  values  C  =  1.001, 1.050, 1.100. 

The  effect  of  the  modified  neutral  equation  of  state  Eq.  (4.1)  on  the  shock  paths  for  the  values  of  a  used 
previously  are  shown  in  Fig.  (4.3)  along  with  the  result  for  C  =  0.  The  apparent  increased  roughness  of  the 
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Fig.  4.2.  shock  path  for  neutral  shock,  shock  with  enhanced  ion  pressure  gradient,  and  shock  with  enhanced  neutral  pressure 
gradient. 


shock  paths  in  Fig.  (4.3)  for  C  =  0  is  an  effect  of  the  y- axis  scale,  which  is  much  finer  in  this  graph  than 
the  others  because  the  shock  is  almost  unperturbed  by  coupling  to  the  charged  particles  alone. 

In  these  examples,  the  pressure  gradient  is  only  enhanced  in  the  charge  separation  region.  The  break 
in  the  curves  when  the  shock  propagation  speed  increase  is  large  is  explained  by  a  feature  of  the  numerical 
method  noted  earlier:  if  the  shock  moves  too  quickly,  the  boundary  conditions  become  contaminated  by 
shock  quantities,  and  the  numerical  solution  is  unable  to  converge. 

The  figures  exhibit  no  increase  of  shock  propagation  speed  when  a  =  10~5,  consistent  with  the  absence  of 
a  charge  separation  region  in  these  cases.  In  fact,  a  small  decrease  in  shock  propagation  speed  is  consistently 
found,  but  this  decrease  is  too  small  to  be  significant.  The  figures  all  show  that  even  if  the  neutral  pressure 
enhancement  is  localized  in  a  small  region,  the  shock  propagation  speed  can  be  increased.  Since  the  pressure 
increase  has  been  chosen  independently  of  the  density  a,  the  increase  in  energy  when  a  =  10”8  is  very  large. 
A  more  useful  way  to  present  the  results  is  to  plot  the  shock  speed  increase  as  a  function  of  some  measure 
of  the  energy  addition  per  particle.  One  possible  measure  is 

C  f 

(4.2)  c  =  - —  /  dxmax{rii(x)  -  ne(x),0} 

lCOLj 


where  the  length  Zc,  defined  by 
(4.3) 


lc 


L 


m{x)>ne{x) 


dx 


is  the  width  of  the  charge  separation  region.  As  noted  above,  lc  ~  Id  «  The  results  are  shown  in  Fig.  (4.4) 
where  the  ‘added  pressure  factor’  is  the  quantity  c  defined  in  Eq.  (4.2). 


5.  Discussion.  Given  that  the  shock  propagation  speed  can  be  increased  if  the  neutral  pressure  can 
be  increased,  even  if  only  in  the  charge  separation  region,  it  is  necessary  to  ask  how  the  neutral  equation 
of  state  could  be  modified  by  non-ideal  gas  effects  on  the  ions;  that  is,  how  can  the  electrostatic  potential 
energy  in  the  charge  separation  region  be  transferred  to  the  neutral  particles? 

Consider  the  model  system  of  an  Enskog  non-ideal  gas  collisionally  coupled  to  an  ideal  gas.  The  equation 
of  state  for  this  system  has  been  computed  by  Thorne  ([11]  Sect.  16.6).  The  pressure  in  the  ideal  gas  is 
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(5.1)  P  =  knT{  1  +  -7i -n'a3x} 
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time  time 


Fig.  4.3.  shock  paths  for  a  =  10“5 [5],  10“ 6 [6],  and  10“ 8 [8]  with  ( upper  left)  C  =  0;  (upper  right)  C  =  1.001;  (lower  left) 
C  —  1.050;  (lower  right)  C  =  1.100.  The  shock  is  accelerated  as  the  strength  of  the  neutral  pressure  gradient  enhancement  C 
increases,  although  the  pressure  enhancement  is  confined  to  the  charge  separation  region  alone. 

where  n/  is  the  number  density  of  Enskog  particles  and  a  is  their  radius,  so  that  n'cr3  is  the  volume  fraction 
occupied  by  the  Enskog  particles.  The  factor  x  is  given  by 

(5.2)  X  =  1  +  J5n>cr3  +  -  -  ’ 

Eq.  (5.1)  shows  that  the  pressure  in  the  ideal  gas  is  enhanced  by  the  presence  of  the  non-ideal  gas.  Since 
long-range  forces  like  the  Coulomb  force  correspond  to  large  a ,  the  effect  of  the  coupling  could  be  large,  even 
if  the  Enskog  gas  is  relatively  rarified. 

Although  this  result  is  suggestive,  its  applicability  to  weakly  ionized  gases  is  doubtful  because  the 
electrostatic  force  does  not  act  on  the  neutral  particles.  To  make  the  model  system  apply  to  a  weakly 
ionized  gas,  the  radius  a  should  take  a  nonzero  value  only  for  ion-ion  or  electron-electron  collisions,  but 
in  ion-neutral  collisions,  a  must  be  set  to  zero.  In  this  case,  the  increased  pressure  predicted  in  Eq.  (5.1) 
disappears. 

Higher  order  interactions,  like  ion-ion-neutral  collisions  will  certainly  be  modified  by  the  presence  of 
the  electrostatic  forces,  but  these  will  not  significantly  change  the  neutral  equation  of  state.  To  conclude, 
although  there  is  considerable  electrostatic  energy  stored  in  the  charge  separation  region,  there  is  no  obvious 
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Fig.  4.4.  shock  propagation  speed  increase  as  a  function  of  energy  added  per  particle 


collisional  mechanism  which  will  transfer  this  energy  to  the  neutral  particles. 

6.  Conclusion.  The  increased  ion  pressure  due  to  non-ideal  gas  effects  in  the  charge  separation  region 
ahead  of  the  shock  front  in  a  weakly  ionized  gas  will  not  significantly  enhance  the  neutral  shock  propagation 
speed.  Other  non-thermal  mechanisms  which  might  lead  to  shock  propagation  speed  enhancement  include 
production  of  molecules  in  excited  states  [12]  and  the  polarizability  of  the  neutral  species  [2],  [13]. 
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